function [img0, img1] = histogram_match(image0, image1);
M = zeros(256,1,'uint8'); %// Store mapping - Cast to uint8 to respect data type
hist1 = imhist(image0); %// Compute histograms
hist2 = imhist(image1);
count_no_zero = sum(hist1>0.5);
display(count_no_zero);
cdf1 = cumsum(hist1) / numel(image0); %// Compute CDFs
cdf2 = cumsum(hist2) / numel(image1);

%// Compute the mapping
for idx = 1 : 256
    [~,ind] = min(abs(cdf1(idx) - cdf2));
    M(idx) = ind-1;
end

%// Now apply the mapping to get first image to make
%// the image look like the distribution of the second image

img0 = M(double(image0)+1);
img1 = image1;

end